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Abstract 

Mixed integer predictive control deals with optimizing integer and real control variables over a receding 
horizon. The mixed integer nature of controls might be a cause of intractability for instances of larger 
dimensions. To tackle this little issue, we propose a decomposition method which turns the original n- 
dimensional problem into n indipendent scalar problems of lot sizing form. Each scalar problem is then 
reformulated as a shortest path one and solved through linear programming over a receding horizon. This 
last reformulation step mirrors a standard procedure in mixed integer programming. The approximation 
introduced by the decomposition can be lowered if we operate in accordance with the predictive control 
technique: i) optimize controls over the horizon ii) apply the first control iii) provide measurement updates 
of other states and repeat the procedure. 



1 Introduction 



Mixed integer predictive control arises when optimizing integer and real control variables in a receding 
horizon context [T]. For this reason, many authors see it as a specific field in the broader area of optimal 
hybrid control [3]. Optimal integer control problems have been receiving a growing attention and are often 
categorized under different names. See, for instance, the literature on finite alphabet control [SJH]. Integer 
control requires a bit more than standard convex optimization techniques. From the literature we know 
that new properties come into play. As an example, look at multimodularity presented as the counterpart 
of convexity in discrete action spaces [3]. When talking about mixed integer variables, it is, of course, not 
possible not to mention the more than vast literature on mixed integer programming [7] . It is exactly in this 
context that we have found inspiration as clarified in more details next. 

In this paper, we have moved our steps along the line of [5] which surveys solution methods for mixed 
integer lot sizing models. Indeed, decomposing an n-dimensional dynamic system into n indipendent lot sizing 
systems is almost all about this paper is centered around. The approximation introduced by the decomposition 
can be reduced if wc operate in accordance with the predictive control technique: i) optimize controls for 
each indipendent system all over a prediction horizon, ii) apply the first control to each indipendent system, 
iii) provide measurement updates of other states and repeat the procedure. The main contribution of this 
work is to reformulate the mixed integer problem of point i) as a shortest path problem and solve this last 
through linear programming. This approach mirrors the method surveyed in with the differences that here 
the shortest path problems run iteratively forward in time over a receding horizon. Reframing the method 
in a receding horizon context is an element of novelty and presents some additional and new issues which are 
discussed and overcome throughout the paper. 
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This paper differs from [T] as we focus on a smaller class of problems that can be solved exactly and do 
not require advanced relaxation methods which, in turn, are a main topic in [1]. To bring our discussion 
back to hybrid control, the lot sizing like model used here has much to do with the inventory example briefly 
mentioned in [3] . There, the authors simply include the example in the large list of hybrid optimal control 
problems but do not address the issue of how to fit general methods to this specific problem. On the contrary, 
this work cannot emphasize enough the computational benefits deriving from the "nice structure" of the lot 
sizing constraints matrix. Binary variables, used to model impulses, match linear programming in a previous 
work of the same author [5]. There, the linear reformulation is a straightforward derivation of the (inverse) 
dwell time conditions appeared first in [6]. Analogies with [2] arc, for instance, the use of total unimodularity 
to prove the exactness of the linear programming reformulation. Differences are in the procedure itself upon 
which the linear program is built up. The shortest path model is an additional element which distinguishes 
the present approach from [5]. 

This paper is organized as follows. We state the problem in Section [21 We then move to present the 
decomposition method in Section [3l In Section |4l we turn to introducing the shortest path reformulation and 
the linear program. We dedicate the last Section [5] to support our theoretical analysis with some numerical 
results. 



2 Mixed integer predictive control 

In mixed integer control we usually have continuous state x{k) G M", continuous controls u{k) G R" and 
disturbances w{k) G M", discrete controls y{k) g {0, 1}" (see e.g., [1]). Evolution of the state over a finite 
horizon of length N is described by a linear discrete time dynamics in the general form ([1]), where A and E 
are matrices of compatible dimensions: 



The above dynamics is characterized by one discrete and continuous control variable per each state, and this 
refiects the idea that we may wish to control indipendently each state component. Also, starting from initial 
state at zero, we wish to drive the final state to zero which is a typical requirement when controlling a system 
over a finite horizon. On this purpose, we have added equality constraints on the final states. Also, we force 
the states to remain confined within a desired region, take for it the positive orthant, which may describe a 
safety region in engineering applications or the desire of preventing shortcomings in inventory applications. 

Continuous and discrete controls are linked together by general capacity constraints ([5]), where the pa- 
rameter C is an upper bound on control: 



For clarity reasons, y{k) is the decision of controUing or not the system, and u{k) is the control action. So if 
we decide not to control the system then the control action is null, otherwise this last is any value between 
zero and its upper bound C. 

The following assumption helps us to describe the common situation where the disturbance seeks to push 
the state out of the desired region. 

Assumption 1 (Unstabilizing disturbance effects) 



x{k + 1) = Ax{k) + Ew{k) + u{k) > 0, a;(0) = x{N) = 0. 



(1) 



< li(fc) < Cy(fc), y(fc)e{0,l}". 



(2) 



Ew{k) < 0. 



(3) 



At this point, the non negative nature of controls u(fc) should become much clearer. Actually, control 
actions are used to push the state far from boundaries into the positive orthant thus to counterbalance the 



unstabilizing effects of disturbances over a certain period to come. However, controlling the system has a 
cost and "over acting" on it is punished by introducing a cost/objective function as explained next. 



The objective function to minimize with respect to y{k) and u{k) is a linear one including proportional, 
holding and fixed cost terms expressed by parameters p'^, h^, and f'^ respectively: 



N-l 



{P'uik) + h^xik) + fy{k)) . (4) 



k=0 



Conditions (II])-(|4l) introduced so far describe coincisely the problem of interest. In the next section, we 
recall a standard method to convert the problem of interest (Il])-(l4|) into a mixed integer linear program 
returning the exact solution in terms of optimal control actions u{k) and y{k). 



Remark 1 For sake of simplicity disturbances w{k) are deterministic and apriori known. The approach 
presented below is still valid if we drop this assumption and turn to consider unknown disturbances. Only, we 
should carefully repropose problem (Qp-Q) in a receding horizon form with iterative measuments updates and 
control optimization forward in time all over the horizon. 



2.1 Mixed integer linear program and exact solution. 



The mixed integer nature of the above program makes it intractable for increasing number of variables and 
horizon length. So, the topic presented below is motivated mainly by comparisons reasons and applies only 
to problems of relatively small dimensions. 



Before introducing the mixed integer linear program we need to define the following notation. Let us start 
by collecting states, continuous and discrete controls, proportional, holding and fixed costs all in opportune 
vectors as shown below: 
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Notice that once we take for and ^/ the value zero, the first and last rows in the aforementioned matrices 
restate the constraints on initial and final state of 



Finally, we are in the condition to establish that problem (IlJ-ljl]) can be solved exactly through the 
following mixed integer linear program: 

{MI PC) min J{u,y) = pu + hx + fy (5) 

Ax + Bm = b (6) 

0<u<Cy, ye {0,1}"^. (7) 



The mixed integer linear program ©-([T]) is the most natural mathematical programming representation 
of the problem of interest (IT])-(|3|). For this reason, throughout this paper we will almost always refer to ©-([7]) 
when we wish to bring back the discussion to the source problem (Il])-(|4]) and its exact solution. 



To overcome the intractability of the mixed integer linear program ©-([T]), we propose a new method 
whose underlying idea is to bring back dynamics ([1} to the lot sizing model [5]. To do this, we introduce 
some additional assumptions on the structure of matrix A which simplify the tractability and affect in no 
way the generality of the results. This argument is dealt with in details in the next section. 



2.2 Introducing some structure on A 

Our main goal in this section is to rewrite ([l} in a "nice" form. With "nice form" we mean a form that 
emphasizes the analogies with standard lot sizing models [5]. "Stop beating around the bush", we will 
henceforth refer to the following dynamics in state of ([T]): 

x{k + 1) = x{k) + Ax{k) + Ew{k) + u(k) > 0. (8) 

The reasons why expression ([5]) is a nice one is that it isolates the dependence of one component state on 
the other ones. To tell it differently we have separated the influence of all other states on state i. It will be 
soon clearer that turning our attention to the new expression ([8]) is a prelude in view of the decomposition 
approach discussed later on. 

Once clarified the reasons, we need next to clarify how to go from ([T]) to ([8]) and what is the underlying 
assumption that allows us to do that. Before doing this let us denote with / G M"^" the identity matrix and 
aij the dependence of state i on state j. So, we can make the following assumption. 



Assumption 2 Matrix A can he decomposed as 



A = / + A, A = 



ai2 . . . ai „_i air, 

021 ... 02 

O-nl 0,n2 ■ ■ ■ 0,n,7i~l 



The reader may notice that ([5]) is a straighforward derivation of ^ once we take for good Assumption O 

Our secondary goal in this section is to preserve the nature of the game which has stabilizing control 
actions playing against unstabilizing disturbances. To do this, in our next assumption we do consider the 
case where the influence of other states on state i is relatively "weak" in comparison to the unstabilizing 
effects of disturbances. 

Assumption 3 (Weakly coupling) 

l^x{k) + Ew{k) <Q. (9) 

Notice that the above assumption preserves the nature of the game by bounding the effects of mutual 
dependence of state components represented by the term Ax{k). A closer look at ([3]) and ([9|) sounds like the 
term Ax{k) do not counterbalance the effects of Ew{k). States mutual dependence only emphasize or reduce 
"weakly" the unstabilizing effects of disturbances. 

We end this section by noticing that ([8|) is not yet in "lot sizing" form [8] . In the next section, we present 
a decomposition approach that translate dynamics ([5]) into n scalar dynamics in "lot sizing" form [5] . 



96 3 Robust decomposition 
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With the term "decomposition" we mean a mathematical manipulation through which the original dynamics 
([5]) is replaced by n independent dynamics of the form: 

x^{k + 1) = x,{k) - d,{k) + u,{k). (10) 

The above dynamics is in a typical lot sizing form in the sense that the (inventory) state tomorrow Xi{k + 1) 
is equal to the (inventory) state today Xi{k) plus the discrepancy between today demand di{k) and today 
reordered quantity Ui{k). Changing ([8]) with (fTO|) is possible once we relate the demand di{k) to the current 
values of all other state components and disturbances as expressed below: 



d^ik) 



E"=i,j5^i A^JX^{k) + ^"^1 E,jWj{k) 
A„xik) + E„w{k)] . 



(11) 



To tell it differently, we do assume that the influence that all other states have on state i enters into equation 
pro through demand di{k) defined in PT|) . Our next step is to make the n dynamics in the form PH)) 
mutually independent. This is possible by replacing the current state values Xj{k), j ^ i with their estimated 
values on the part of agent i which we denote by Xjik), j ^ i. Still with reference to (|10p. this implies to 
replace the current demand di(k) by the "estimated" demand d,;(fc) defined as in ((T2|) where X'' is the set of 
admissible state vectors x{k): 

d^{k) = max {-A,.e - E,,w{k)) . (12) 

The idea behind (fT^ is to take for estimated value the worst admissible demand, i.e., the demand that would 
push the state out of the positive orthant in a fewest time and such a demand is of course the maximal one. 
However, it must be noted that we cannot see any drawbacks in combining other decomposition methods 
with the approach presented in the rest of the paper. To complete the decomposition, it is left to turn the 
objective function (j4]) into n indipcndent components 

N-l 



J,iu„ y,) = (pN'(^) + h^,x^{k) + f^y^{k)) 



k=0 

Note that because of the linear structure of J{u,y) in ([5]), it turns J{u,y) = ^"^i Jiiui,yi). So, in the 
end we have translated our original problem into n indipcndent mixed integer linear minimization problems 
of the form p3p - (fT5l) as requested at the beginning of this section. In the spirit of predictive control, each 
minimization problem is then solved forwardly in time all over the horizon. So, for r = 0, . . . , iV — 1 we need 
to solve 



(M/PC) min V (pS^(^) + /ifx,(A:) + /fy,(fc)) (13) 

X.,{k + 1) ^ X,{k) ^ d^{k) + U^{k) > 0, x,{t)=^°,X,{N) = (14) 

0<u,{k) <Cy,ik), y,{k)e {0,1}. (15) 

It is worth to be noted that non null initial states, which materialize in values of strictly greater than 
zero in constraints p4|l might induce infcasibility of (MlPd). So. moving from {MIPC) to {MlPd) has 
this little drawback that we will discuss in more details later on in Section together with some other issues 
concerned with the receding implementation of our method. 



4 Shortest path and linear programming 

So far, we have first formulated the problem of interest and then decomposed it into n indipendent scalar 
problems. By the way, decomposition is only the first step of our solution approach. Actually, the mixed 



116 integer nature of variables in p^ - p^ is still an issue to be dealt with. This second part of the work focuses on 

117 the relaxation of the integer constraints yi{k) G {0, 1} which would facilitate the tractability of the problem. 

118 It is well known that relaxation introduces, in general, some approximation in the solution. The main result 

119 of this work establishes that, for the problem at hand, relaxing and massaging the problem in a certain 

120 manner, will lead to a shortest path reformulation of the original problem. This is a great result as, it is 

121 well known that shortest path problem are in turn easily tractable and solvable through linear programming. 

122 Shortest path formulations are based on the notion of regeneration interval discussed in details in the next 

123 section. 



124 4.1 Regeneration interval [a, (3] 



125 Let us start by introducing a formal definition of regeneration interval which represents the central topic in 

126 this section. The definition, available in the literature for scalar lot sizing models, is borrowed from [5] and 

127 adapted to each single (scalar) dynamics i of our decomposed n-dimensional model. So, with reference to 

128 the generic minimization problem i expressed by (|13p - (|15p . let us state what follows. 



129 Definition 1 (Pochet and Wolsey 1993) A pair of periods [a, [3] form a legeneiationinteival for {xi,Ui,yi) 

130 if Xi{a 1) = Xi{(3) = and Xi{k) > for fc = a, a + 1, . . . , /3 — 1. 



131 Given a regeneration interval [a,/3], we can define the accumulated demand over the interval df , and 

132 the residual demand rf ^ as 



^d~(fc). 



aP 



jaP 



c 



c. 



(16) 



133 Our idea is now to translate problem ((T3| - (|15p into new variables. More formally, let us consider variables 

134 yf^{k) and e"^{k) defined in (|17p with the following meaning. Variable yf^{k) is equal to one in presence of 

135 a saturated control on time k and zero otherwise. Similarly, variable e"^(fc) is equal to one in presence of a 

136 non saturated control on time k and zero otherwise: 

^ f 1 if u,{k) = C ^ap,^. ^ f 1 if < u,{k) < C 
' ^ [0 otherwise. i \ f |^ q otherwise. 

137 To translate the meaning of yf^ {k) and e"'^(fc) in a lot sizing context, such variables tell us on which period 

138 full or partial batches are ordered. 

At this point and with in mind the above variable transformation, we can rely on well known results in 
the lot sizing literature which convert the original mixed integer problem (|13 p -(|15 p into a number of linear 

programs (^LPt^) , each one associated to a specific regeneration interval. Regeneration intervals and the 
associated linear programs are mutually related in a way that gives raise to a shortest path problem, which 
will be the central topic in the next section. For now, we simply reproposc below the linear programming 
problem associated to a single regeneration interval [a, /?]. Denoting by e'^ = pf + J^j'^k+i ^""^"^ after some 



standard manipulation, the linear program for fixed regeneration interval [a, /3] appears as: 



^LP"^^ mm 
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t = a, ...,/? -f 



yfik),efik)>0, k^a,...,p. 



(18) 
(19) 
(20) 

(21) 

(22) 
(23) 



The above model is extensively used in the lot sizing context. We can limit ourselves to a pair of comments 
on the underlying idea of the constraints. So, let us start by focusing on the equality constraints (|19p and 
(PT|) . These constraints tell us that the ordered quantity over the interval has to be equal to the accumulated 
demand over the same interval. This makes sense as initial and final state of a regeneration interval are null 
by definition. Let us turn our attention to the inequality constraints (|20p and (|22p . There, we impose that 
the accumulated demand in any subinterval may not exceed the ordered quantity over the same subinterval. 
Again, this is due to the condition that states are nonnegative at any period of a regeneration interval. 
Finally, the objective function (|18l) is simply a rearrangement of (|13p induced by the variable transformation 
seen above and specialized to the regeneration interval [a, /3] rather than on the entire horizon [0, N]. 



We are ready to recall the following "nice property" of {LP"^) presented first by Pochet and Wolsey in 



150 Theorem 1 (Total unimodularity) The optimal solution of {LP"^) is feasible 



Proof. The proof is based on the observation that the constraint matrix of {LP°'^) is a — 1 matrix. We 
can reorder the constraints in a certain manner, so that matrix has the consecutive I's property on each 
column and turns to be totally unimodular. It follows that y"'^ and e"'^ are — 1 in any extreme solution. 



The above theorem represents a first step in the process of converting the mixed integer problem [MlPCi] 
into a linear programming one. 



157 4.2 Shortest path 



In the previous section we have introduced a linear programming problem associated to a specific regeneration 
interval. In this section, we resort to well known results on lot sizing to come up with a shortest path model 
which links together the linear programming problems of all possible regeneration intervals. Actually, it must 
be noted that the solution of (|13p - (fTS)) can be expressed as a unique regeneration interval [0, A^] or as a list 
of regeneration intervals. 



163 So, let us define variables z,^ 

164 appears or not in the solution of 

165 form below. For t = 0, . . . , A^ — 1 



"'^ G {0, 1} which tell us one or zero whenever a regeneration interval [a, (3] 
O]) - p3|) . The linear programming problem solving -(US]) takes on the 
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Let us spend a couple of words on the meaning of the above hnear program. Constraints (|27p - (PT|) should 
be familiar to the reader as they already appeared in ([T^ - (P5)) . The only difference is that, now, because of 
the presence of zf in the right hand term, the constraints referring to a given regeneration interval come 
into play only if that interval is chosen as part of the solution, that is, whenever zf is set equal to one. 
Furthermore, a new class of constraints appear in (|25 |) - (f26| . These constraints are typical of shortest path 
problems and in this specific case help us to force the variables zf{k) to describe a path from to N. 
Finally, note that for r = 0, the linear program (LPi) coincide with the linear program presented by Pochct 
and Wolsey in . 



174 At this point, we are in a position to recall the crucial result established in [8]. 



175 Theorem 2 (Pochet and Wolsey, 1993) The linear program (LP^) solves [MlPCi). 



176 Proof. (Sketch) It turns out that the linear program (LP^) is a shortest path problem on variables . Arcs 

177 are all associated to a different regeneration interval [a, /?] and the respective costs are the optimal values of 

178 the objective functions of the corresponding linear programs (LP"'^). We refer the reader to [5] for further 

179 details. 

□ 



180 4.3 Receding horizon implementation of (LPj) 



This section is dedicated to certain issues concerned with the implementation of iLPi) in a receding horizon 
context as typical of predictive control. As the reader may know, in predictive control we solve iLPi) 
iteratively and forward in time all over the horizon. In the formulation of (LPi), this is stated clearly when 



we specify that r goes from to — 1 and for each value of r we obtain a new linear program of type (LPi). 
After we solve (LPi) for t = 0, we apply the first control to the system, update initial states according to 
the last available measurements at time t ~ I and move to solve a new (LPi) starting at r = 1. We repeat 
this procedure until the end of the horizon, t = TV — 1. So, consecutive linear programs are linked together 
by initial state condition expressed in ([T4| . and which we rewrite below 

181 At this point, we would restate with emphasis the fact that dealing with non null initial states is a main 

182 difference between the linear program (LPi) and the linear program used in the lot sizing literature [5]. To 

183 counter this little issue, we need to elaborate more on how to compute the accumulated demand in (|16p . 

184 Actually, take for [r, i] any interval with x(r) = > sO. Then, condition ([TC]) needs to be revised as 

<*=™^^|E^~'W-^?'o|- (32) 

185 The rational behind the above formula has an immediate interpretation in the lot sizing context. Actually, 

186 the effective demand over an interval is the accumulated demand reduced by the inventory stored and initially 

187 available at the warehouse. From a computational standpoint, the revised formula ([5^ has a different effect 

188 depending on the cases where the accumulated demand exceeds the initial state or not as discussed next. 



1- Sfe=Q d.i{k) > the mixed linear program [MPCi) with initial state x[t) = > and accumulated 

190 demand Yl,k=a di{k) is turned into an (LPi) characterized by null initial state x{a — 1) = and effective 

191 demand d"^ — X]fc=a d,i{k) — ^.^ as in the example below: 

(MPCi) ^d,(fc) = 12, x(r)==$° = 10 =^ {LPi) x{a-l) = 0, df = 2; 

k—a 

2. J2k=a d-i{k) < the mixed linear program (MPd) with initial state x{t) = > and accumulated 

193 demand X]fc=ct^j(^) unfeasible. The solution obtained at previous period r — 1 applies. A second 

194 example is shown next: 

/3 

(MPC,) ^di{k) = 7, a;(r)=^° = 10 =^ (LP,) unfeasible. 

k—a 

195 In both cases, the revised formula ((32|) helps us to generalize the linear program (LPi) to cases where 
19b the initial state is non null and this is a crucial point when applying the lot sizing model in a receding 
197 horizon form. 



198 5 Numerical example 



In this specific example, dynamics ([IJ takes on the form expressed below. Such a dynamics is particularly 
significative as it reproduces the typical influence between position and velocity in a sampled second-order 
system. Initial and final states are null and state values must remain in the positive quadrant all over the 
horizon. More specifically, denoting by xi the position and X2{k) an opposite in sign velocity, the dynamics 
appears as: 



xi [k 

X2{k 



K 1 



xi{k) 

X2{k) 



wi{k) 

W2{k) 



ui(k) 

U2{k) 



>0, 



a:i(0) 

X2(0) 



XI (TV) 

X2{N) 



0. (33) 



A closer look at the first equation reveals that a greater velocity X2{k) reflects into a faster decrease of position 
xi{k + 1). Similarly, the second equation tells us that a greater position xi{k) induces a faster increase of 



206 velocity 0:2(^+1) because of some elastic reaction. In both equations, the non negative disturbances Wi{k) < 

207 seek to push the states a;,;(fc) out of the positive quadrant in accordance to Assumption [31 Their effect is 

208 counterbalanced by positive control actions Ui. Notice that matrix A can be decomposed as described 

209 in Assumption [2] Also, acting on parameter k, we can easily guarantee the "weakly coupling" condition 

210 expressed in Assumption [3l 



Turning to the capacity constraints for this two-dimensional example, these constraints can be rewrit- 
ten as: 

U2(k) J [ y2[k) J [ y2{k) J 

It is left to comment on the objective function (|4]). We consider the case where fixed costs are much more 
relevant than proportional and holding ones. This materializes in choosing a high value for f'' in comparison 
to values of parameters p'^, as shown in the next linear objective function: 

N-l 



A:=0 



This choice makes sense for two reasons. First, all the work is centered around issues deriving from the 
integer nature of y{k). So, high values of emphasize the role of integer variables in the objective function. 
Second, high fixed costs incentivate solutions with the fewest number of control actions and this facilitate 
the validation and interpretation of the simulated results. 



The next step is to decompose dynamics (|33p in scalar lot sizing form (|14p which we rewrite below: 

Xi{k + 1) = Xi{k) — di{k) + Ui{k). 

215 When it comes to the discussion on how to compute the estimated demand d^, a natural choice is to set di as 

216 below, where we have denoted by ii(fc) (respectively X2{k)) the estimated value of state xi{k) (respectively 

217 X2(k)) available to agent 2 (agent 1): 
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(34) 



218 Now, the question is: which expression should we use to represent the set of admissible state vectors 

219 appearing in equation p^ ? This question has much to do with another one: how does agent 1 predict X2 

220 and the same for agent 2 with respect to state xi? A possible answer is shown next: 
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(35) 



221 Let us elaborate more on the above equations. Regarding to variable X2(fc), this is used in the evolution of 

222 di{k) as in the first equation of (p4l) . Because of the positive contribution of the term KX2{k) on di{k), a 

223 conservative approach would suggest to take for X2{k) a possible upper bound of X2{k) and this is exactly 

224 the spirit behind the evolution of 2:2 (fc) as expressed in the second equation of psp . Here, xi is an average 

225 value for xi. A similar reasoning applies to ii(fc), used in the evolution of d2(k) as in the second equation of 

226 (p4l) . Wc now observe a negative contribution of the term — Kii(fc) on d2{k) and therefore take for xi(k) a 

227 possible lower bound of xi{k) as shown in the first equation of ([35]). 



228 We can now move to show and comment our simulated results. We have carried out two different set 

229 of experiments whose parameters are displayed in Table [TJ In the line of the weakly coupling assumption 

230 (sec Assumption [3|) , wc have set k, small enough and in the range equal from 0.01 to 0.225. Such a range 

231 works good as wc will sec that \K.Xi\ is always less than Wi, which also means Ax'(fc) -f- Ew{k) < 0. For sake 

232 of simplicity and without loss of generality, capacity C is set to three, disturbances Wi arc unitary and xi 

233 is equal to one. Unitary disturbances facilitate the check out and interpretation of the results as when the 

234 accumulated demand over the horizon turns to be very close to the horizon length. The two experiments 

235 differ also in the horizon length N for the reasons clarified next. 
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K 
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Wi(fc) 


W2{k) 




I 


1 ...10 


0.1 


3 


1 


1 


1 


II 


6 


{0.01, 0.2, 0.225} 


3 


1 


1 


1 



Table 1: Simulation parameters chosen for the two experiments. 



236 The first set of experiments aims at analysing the computational benefits of decomposition and relaxation 

237 upon which our solution method is based. So, we consider horizon lenghts N from one to ten. We do not need 

238 to consider larger values of N as even in this small range of values, differences in the computational times 

239 are already evident enough as clearly illustrated in Fig. [T] Here, we plot the average computational time vs. 

240 the horizon length N of the mixed integer predictive control problem (solid diamonds), of the decomposed 

241 problem (MlPCi) (dashed squares), and of the linear program (LPi). Average computational time means 

242 the average time for one agent to make a single decision (the total time is about 2A'' times the average one) . 

243 As the reader may notice, the computational time of the linear program (LPi) is a fraction either of the one 

244 requested by the (MPC) or of the one required by the [MlPCi). 

245 (Figure [1] about here) 

In a second set of simulations, wc have inspected how the percentage error 

^ optimal cost of (MPd) — optimal cost of (MPC) ^ 
optimal cost of (MPC) 

246 varies with different values of the elastic coefficient k. The role of k is crucial as we recall that k describes the 

247 effective tightness and coupling between different states xi{k) and X2{k). We do expect that small values for 

248 coefficient k, which means weak coupling of state components, may lead to small errors e%. Differently, high 

249 values of k, describing a strong coupling between state components, are supposed to induce higher values of 

250 e%. 

251 This is in line with what we can observe in Fig. |3] where we plot the error e% as function of coefficient k. 

252 For a relatively small values of k in the range from to 0.2, we observe a percentage error not exceeding the 

253 one percent, e% < 1. A discountinuity at around k = 0.2 causes the error e% to go from about 1% to 20%. 

254 (Figure [2] about here) 

255 We might not be surprised as discountinuity of errors is typical in mixed integer programs and wc try to 

256 clarify this in more details in the plot of Fig. 31 Here, for a horizon length = 6 and for a relatively high 

257 value of K = 0.225; we display the exact solution (dashed squares) and approximate solution (solid triangles) 

258 returned by the mixed integer linear program (MIPC) and by the linear program (LPi) respectively. The 

259 solution is in terms of the time plot of states Xi{k), continuous controls Ui{k) and discrete controls yi{k). 

260 Dotted lines represent predicted trajectories in earlier periods of the receding horizon implementation. At a 

261 first check, and this is in accordance with what we do expect, we note that controls Ui(k) never exceed the 

262 capacity and are always associated to unitary control actions yi{k). Now, with a look at the behaviour of 

263 discrete controls yi{k), it can be observed that the approximate solution presents four control actions (four 

264 peaks at one), whereas the exact solution has control yi{k) acting on the system only three times (three 

265 peaks at one). One peak out of four represents an increase in the use of control actions of about 25 percent 

266 which reflects into an approximate increase in the percentage error of 20%. A last observation concerning 

267 the exact plot of yi{k) is that the number of control actions are as minimal as possible, i.e., three for yi{k) 

268 and two for y2{k). This makes sense as the accumulated demand over the horizon approximates by above the 

269 horizon length. This implies that the minimum number of control actions can be roughly obtained dividing 

270 the accumulated demand (about something above six) by the capacity C (equal to three) and rounding the 

271 fractional result up to the next integer. 



272 (Figure 13] about here) 

273 Let us move to compare exact and approximate solutions for a smaller value of k = 0.2. With reference 

274 to Fig. 01 we observe that, differently from above, discrete controls yi{k) coincide. However, we still have 

275 notable differences in the plot of continuous controls ui{k) which cause distinct state trajectories for Xi{k). 

276 Small differences can be noted for U2{k) and X2{k) as well. The observed differences still cause a reduced 

277 percentage error e% = 1. 

278 (Figure H] about here) 

279 We conclude our simulations by showing that the percentage error e% is around zero when we reduce further 

280 the value of k to 0.01. This is evident if we look at Fig. \E\ where plots of different styles overlap which means 

281 that exact and approximate solutions coincide. 

282 (Figure [5] about here) 
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Figure 2: Percentage error e% for different values of the elastic coefficient k. 
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Figure 3: Elastic coefHcient k = 0.225. Exact solution (dashed squares) and approximate solution (solid 
triangles) returned by the mixed integer linear program (MI PC) and by the linear program [LPi) respectively. 
Horizon length TV = 6. Time plot of states Xi(k), continuous controls Ui{k) and discrete controls yi{k). 
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Figure 4: Elastic coefRcient k, = 0.20. Exact solution (dashed squares) and approximate solution (solid 
triangles) returned by the mixed integer linear program [MI PC) and by the linear program {LPi) respectively. 
Horizon length iV = 6. Time plot of states Xi{k), continuous controls Ui{k) and discrete controls yi{k). 
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Figure 5: Elastic coefficient k = 0.001. Exact solution (dashed squares) and approximate solution (solid 
triangles) returned by the mixed integer linear program (MI PC) and by the linear program (LPi) respectively. 
Horizon length iV = 6. Time plot of states Xi{k), continuous controls Ui{k) and discrete controls yi{k). 
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